**MC models**
** this same code can be used, with small modification, for the 30 and 50 percent cutoffs as well.
clear all
set more off
cap prog drop mclogit
cap prog drop mclogit_p
cap prog drop mcre
cap prog drop mcre_p
cap log cl
set seed 2356935

adopath + "${user}Code\"

* Read-in land use data *
******************************************************************************
use "${data}BugandaAnalysisME.dta", clear

la var zero "2008-2009" 

replace ruggedness = ln(1+ruggedness)
replace maize = maize/100
replace elevation = ln(1+elevation)


keep if gte10 == 1

drop if period == 2000
qui {
	sort id year
	xtset id

	*define vars to be used
	loc y "binhan10"	
	loc x0 "pi1 pi2 pi3 pi5 pi6 pi7 pi8 pi9 "
	loc x1 " mailo1 mailo2 mailo3 mailo5  mailo6 mailo7 mailo8 mailo9  free1 free2 free3 free5 free6 free7 free8 free9 maize ruggedness dist_kampa elevation"
	loc x1m "maize ruggedness dist_kampa "

	
	g z2=num_scenes 
	g z1=num_scenes *ruggedness
	lab var z2 "Cloud-free scenes, L7"
	lab var z1 "Cloud-free scenes, L7 x ruggedness"

	loc z1 "z1"
	loc z2 "z2"

	*create means for CRE models*
	loc mx1 " "
	loc mz1 " "
	loc mz2 " "
	foreach x in `x1m' {
		bys code3: egen m`x'=mean(`x')
		loc mx1 = "`mx1' m`x'"
	}

	foreach x in `z2' {
		bys code3: egen m`x'=mean(`x')
		loc mz2 = "`mz2' m`x'"
	}

}

*estimators: misclassification*
*mc logit**
loc i=1

logit `y' `x1' `x0' `mx1', vce(cl code3) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in `x1' `mx1' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}

mclogit `y' `x1' `mx1', vce(cl code3) mc(`z1' ) diff from(a0:_cons=-3 a0:z1=-.01 a1:_cons=1.5 a1:z1=0.001 `fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in `x1' `mx1'  {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}

foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}
mclogit `y' `x0' `x1' `mx1', vce(cl code3) mc(`z1') diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in `x0' `x1' `mx1' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}
foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}

*from(a0:z2 = -0.005 a1:z2=-0.04 `fromp')
mclogit `y' `x0' `x1' `mx1', vce(cl code3) mc(`z1' `z2') diff from( `fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
g G0p=[a0]_cons
g G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)
eststo col`i': margins, dydx(`x1') predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')

loc i=`i'+1

drop G0p G1p 

esttab col*  ///
	using "${tables}BugandaMEResultsLogitgte10.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" ) ///
	 keep( `x1') nonumbers  replace  fragment  ///
	 scalar(G0 G1 logL) prehead({ \begin{tabular}{lc}    ///
	 \hline \hline ///
	 \\  ///
	 &  (2)  \\ ) ///
	 prefoot( \\ ) ///
		postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate sample. ///
		CRE means included at district level.  Standard errors are clustered at the county level. ///
		Marginal effects evaluated at sample means are displayed. $G_{0}$ and $G_{1}$ are the probability of a false positive and negative, respectively, ///
		evaluated at sample means. The false positive/negative rates to depend upon the number of L7 cloud-free scenes ///
		and the interaction between cloud-free scenes and ruggedness. Time fixed effects included in ///
		all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

	